// @file nbtheory.cpp This code provides number theory utilities.
// @author TPOC: contact@palisade-crypto.org
//
// @copyright Copyright (c) 2019, New Jersey Institute of Technology (NJIT)
// All rights reserved.
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
// 1. Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution. THIS SOFTWARE IS
// PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
// EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#define _USE_MATH_DEFINES

#include <time.h>
#include <cmath>

#include <chrono>
#include <sstream>

#include "math/distributiongenerator.h"
#include "math/nbtheory.h"

#include "utils/debug.h"

namespace lbcrypto {

/*
 Generates a random number between 0 and n.
 Input: BigInteger n.
 Output: Randomly generated BigInteger  between 0 and n.
 */
template <typename IntType>
static IntType RNG(const IntType &modulus) {
  // static parameters for the 32-bit unsigned integers used for multiprecision
  // random number generation
  static const usint chunk_min = 0;
  static const usint chunk_width = std::numeric_limits<uint32_t>::digits;
  static const usint chunk_max = std::numeric_limits<uint32_t>::max();

  static std::uniform_int_distribution<uint32_t> distribution(chunk_min,
                                                              chunk_max);

  // Update values that depend on modulus.
  usint modulusWidth = modulus.GetMSB();
  // Get the number of chunks in the modulus
  // 1 is subtracted to make sure the last chunk is fully used by the modulus
  usint chunksPerValue = (modulusWidth - 1) / chunk_width;

  // result is initialized to 0
  IntType result;

  // temp is used for intermediate multiprecision computations
  IntType temp;

  // stores current random number generated by built-in C++ 11 uniform generator
  // (used for 32-bit unsigned integers)
  uint32_t value;

  do {
    result = 0;
    // Generate random uint32_t "limbs" of the BigInteger
    for (usint i = 0; i < chunksPerValue; i++) {
      // Generate an unsigned long integer
      value = distribution(PseudoRandomNumberGenerator::GetPRNG());
      // converts value into IntType
      temp = value;
      // Move it to the appropriate chunk of the big integer
      temp <<= i * chunk_width;
      // Add it to the current big integer storing the result
      result += temp;
    }

    // work with the remainder - after all 32-bit chunks were processed
    temp = modulus >> chunksPerValue * chunk_width;

    // Generate a uniform number for the remainder
    // If not 1, i.e., the modulus is either 1 or a power of 2*CHUNK_WIDTH
    if (temp.GetMSB() != 1) {
      uint32_t bound = temp.ConvertToInt();

      // default generator for the most significant chunk of the multiprecision
      // number
      std::uniform_int_distribution<uint32_t> distribution2 =
          std::uniform_int_distribution<uint32_t>(chunk_min, bound);

      value = distribution2(PseudoRandomNumberGenerator::GetPRNG());
      // converts value into IntType
      temp = value;
      // Move it to the appropriate chunk of the big integer
      temp <<= chunksPerValue * chunk_width;
      // Add it to the current big integer storing the result
      result += temp;
    }
  } while (result > modulus);  // deals with the rare scenario when the bits in
                               // the most significant chunk are the same
  // and the bits in the following chunk of the result are larger than in the
  // modulus
  return result;
}
/*
 A witness function used for the Miller-Rabin Primality test.
 Inputs: a is a randomly generated witness between 2 and p-1,
 p is the number to be tested for primality,
 s and d satisfy p-1 = ((2^s) * d), d is odd.
 Output: true if p is composite,
 false if p is likely prime
 */
template <typename IntType>
static bool WitnessFunction(const IntType &a, const IntType &d, usint s,
                            const IntType &p) {
  DEBUG_FLAG(false);
  DEBUG("calling modexp a " << a << " d " << d << " p " << p);
  IntType mod = a.ModExp(d, p);
  DEBUG("mod " << mod);
  bool prevMod = false;
  for (usint i = 1; i < s + 1; i++) {
    DEBUG("wf " << i);
    if (mod != IntType(1) && mod != p - IntType(1)) {
      prevMod = true;
    } else {
      prevMod = false;
    }
    mod.ModMulFastEq(mod, p);
    if (mod == IntType(1) && prevMod) {
      return true;
    }
  }
  return (mod != IntType(1));
}

/*
 A helper function to RootOfUnity function. This finds a generator for a given
 prime q. Input: BigInteger q which is a prime. Output: A generator of prime q
 */
template <typename IntType>
static IntType FindGenerator(const IntType &q) {
  DEBUG_FLAG(false);
  std::set<IntType> primeFactors;
  DEBUG("FindGenerator(" << q << "),calling PrimeFactorize");

  IntType qm1 = q - IntType(1);
  IntType qm2 = q - IntType(2);
  PrimeFactorize<IntType>(qm1, primeFactors);
  DEBUG("prime factors of " << qm1);
#if !defined(NDEBUG)
  for (auto &v : primeFactors) DEBUG(v << " ");
#endif
  bool generatorFound = false;
  IntType gen;
  while (!generatorFound) {
    usint count = 0;

    // gen = RNG(qm2).ModAdd(IntType::ONE, q); //modadd note needed
    gen = RNG(qm2) + IntType(1);
    DEBUG("generator " << gen);
    DEBUG("cycling thru prime factors");

    for (auto it = primeFactors.begin(); it != primeFactors.end(); ++it) {
      DEBUG(qm1 << " / " << *it << " " << gen.ModExp(qm1 / (*it), q));

      if (gen.ModExp(qm1 / (*it), q) == IntType(1))
        break;
      else
        count++;
    }
    if (count == primeFactors.size()) generatorFound = true;
  }
  return gen;
}

/*
 A helper function for arbitrary cyclotomics. This finds a generator for any
 composite q (cyclic group). Input: BigInteger q (cyclic group). Output: A
 generator of q
 */
template <typename IntType>
IntType FindGeneratorCyclic(const IntType &q) {
  DEBUG_FLAG(false);
  std::set<IntType> primeFactors;
  DEBUG("calling PrimeFactorize");

  IntType phi_q = IntType(GetTotient(q.ConvertToInt()));
  IntType phi_q_m1 = IntType(GetTotient(q.ConvertToInt()));

  PrimeFactorize<IntType>(phi_q, primeFactors);
  DEBUG("done");
  bool generatorFound = false;
  IntType gen;
  while (!generatorFound) {
    usint count = 0;
    DEBUG("count " << count);

    gen = RNG(phi_q_m1) + IntType(1);  // gen is random in [1, phi(q)]
    if (GreatestCommonDivisor<IntType>(gen, q) != IntType(1)) {
      // Generator must lie in the group!
      continue;
    }

    // Order of a generator cannot divide any co-factor
    for (auto it = primeFactors.begin(); it != primeFactors.end(); ++it) {
      DEBUG("in set");
      DEBUG("divide " << phi_q << " by " << *it);

      if (gen.ModExp(phi_q / (*it), q) == IntType(1))
        break;
      else
        count++;
    }

    if (count == primeFactors.size()) generatorFound = true;
  }
  return gen;
}

/*
 A helper function for arbitrary cyclotomics. Checks if g is a generator of q
 (supports any cyclic group, not just prime-modulus groups) Input: Candidate
 generator g and modulus q Output: returns true if g is a generator for q
 */
template <typename IntType>
bool IsGenerator(const IntType &g, const IntType &q) {
  DEBUG_FLAG(false);
  std::set<IntType> primeFactors;
  DEBUG("calling PrimeFactorize");

  IntType qm1 = IntType(GetTotient(q.ConvertToInt()));

  PrimeFactorize<IntType>(qm1, primeFactors);
  DEBUG("done");

  usint count = 0;

  for (auto it = primeFactors.begin(); it != primeFactors.end(); ++it) {
    DEBUG("in set");
    DEBUG("divide " << qm1 << " by " << *it);

    if (g.ModExp(qm1 / (*it), q) == IntType(1))
      break;
    else
      count++;
  }

  if (count == primeFactors.size())
    return true;
  else
    return false;
}

/*
 finds roots of unity for given input.  Assumes the the input is a power of two.
 Mostly likely does not give correct results otherwise. input:  m as number
 which is cyclotomic(in format of int), modulo which is used to find generator
 (in format of BigInteger)

 output:  root of unity (in format of BigInteger)
 */
template <typename IntType>
IntType RootOfUnity(usint m, const IntType &modulo) {
  DEBUG_FLAG(false);
  DEBUG("in Root of unity m :" << m << " modulo " << modulo.ToString());
  IntType M(m);
  if ((modulo - IntType(1)).Mod(M) != IntType(0)) {
    std::string errMsg =
        "Please provide a primeModulus(q) and a cyclotomic number(m) "
        "satisfying the condition: (q-1)/m is an integer. The values of "
        "primeModulus = " +
        modulo.ToString() + " and m = " + std::to_string(m) +
        " do not satisfy this condition";
    PALISADE_THROW(math_error, errMsg);
  }
  IntType result;
  DEBUG("calling FindGenerator");
  IntType gen = FindGenerator(modulo);
  DEBUG("gen = " << gen.ToString());

  DEBUG("calling gen.ModExp( "
        << ((modulo - IntType(1)).DividedBy(M)).ToString() << ", modulus "
        << modulo.ToString());
  result = gen.ModExp((modulo - IntType(1)).DividedBy(M), modulo);
  DEBUG("result = " << result.ToString());
  if (result == IntType(1)) {
    DEBUG("LOOP?");
    result = RootOfUnity(m, modulo);
  }

  /*
   * At this point, result contains a primitive root of unity. However,
   * we want to return the minimum root of unity, to avoid different
   * crypto contexts having different roots of unity for the same
   * cyclotomic order and moduli. Therefore, we are going to cycle over
   * all primitive roots of unity and select the smallest one (minRU).
   *
   * To cycle over all primitive roots of unity, we raise the root of
   * unity in result to all the powers that are co-prime to the
   * cyclotomic order. In power-of-two cyclotomics, this will be the
   * set of all odd powers, but here we use a more general routine
   * to support arbitrary cyclotomics.
   *
   */

  IntType mu = modulo.ComputeMu();
  IntType x(1);
  x.ModMulEq(result, modulo, mu);
  IntType minRU(x);
  IntType curPowIdx(1);
  std::vector<IntType> coprimes = GetTotientList<IntType>(m);
  for (usint i = 0; i < coprimes.size(); i++) {
    auto nextPowIdx = coprimes[i];
    IntType diffPow(nextPowIdx - curPowIdx);
    for (IntType j(0); j < diffPow; j += IntType(1)) {
      x.ModMulEq(result, modulo, mu);
    }
    if (x < minRU && x != IntType(1)) {
      minRU = x;
    }
    curPowIdx = nextPowIdx;
  }
  return minRU;
}

template <typename IntType>
std::vector<IntType> RootsOfUnity(usint m, const std::vector<IntType> moduli) {
  std::vector<IntType> rootsOfUnity(moduli.size());
  for (usint i = 0; i < moduli.size(); i++) {
    rootsOfUnity[i] = RootOfUnity(m, moduli[i]);
  }
  return rootsOfUnity;
}

template <typename IntType>
IntType GreatestCommonDivisor(const IntType &a, const IntType &b) {
  DEBUG_FLAG(false);
  IntType m_a, m_b, m_t;
  m_a = a;
  m_b = b;
  DEBUG("GCD a " << a << " b " << b);
  while (m_b != IntType(0)) {
    m_t = m_b;
    DEBUG("GCD m_a.Mod(b) " << m_a << "( " << m_b << ")");
    m_b = m_a % m_b;

    m_a = m_t;
    DEBUG("GCD m_a " << m_b << " m_b " << m_b);
  }
  DEBUG("GCD ret " << m_a);
  return m_a;
}

/*
 The Miller-Rabin Primality Test
 Input: p the number to be tested for primality.
 Output: true if p is prime,
 false if p is not prime
 */
template <typename IntType>
bool MillerRabinPrimalityTest(const IntType &p, const usint niter) {
  DEBUG_FLAG(false);
  if (p < IntType(2) || ((p != IntType(2)) && (p.Mod(2) == IntType(0))))
    return false;
  if (p == IntType(2) || p == IntType(3) || p == IntType(5)) return true;

  IntType d = p - IntType(1);
  usint s = 0;
  DEBUG("start while d " << d);
  while (d.Mod(2) == IntType(0)) {
    d.DividedByEq(2);
    s++;
  }
  DEBUG("end while s " << s);
  bool composite = true;
  for (usint i = 0; i < niter; i++) {
    DEBUG(".1");
    IntType a = RNG(p - IntType(3)).ModAdd(2, p);
    DEBUG(".2");
    composite = (WitnessFunction(a, d, s, p));
    if (composite) break;
  }
  DEBUG("done composite " << composite);
  return (!composite);
}

/*
 The Pollard Rho factorization of a number n.
 Input: n the number to be factorized.
 Output: a factor of n.
 */
template <typename IntType>
const IntType PollardRhoFactorization(const IntType &n) {
  DEBUG_FLAG(false);
  IntType divisor(1);

  IntType c(RNG(n));
  IntType x(RNG(n));
  IntType xx(x);

  // check divisibility by 2
  if (n.Mod(2) == IntType(0)) return IntType(2);

  // Precompute the Barrett mu parameter
  IntType mu = n.ComputeMu();

  do {
    x = x.ModMul(x, n, mu).ModAdd(c, n, mu);
    xx = xx.ModMul(xx, n, mu).ModAdd(c, n, mu);
    xx = xx.ModMul(xx, n, mu).ModAdd(c, n, mu);

    divisor = GreatestCommonDivisor((x > xx) ? x - xx : xx - x, n);

    DEBUG("PRF divisor " << divisor.ToString());
  } while (divisor == IntType(1));

  return divisor;
}

/*
 Recursively factorizes and find the distinct primefactors of a number
 Input: n is the number to be prime factorized,
 primeFactors is a set of prime factors of n.
 */
template <typename IntType>
void PrimeFactorize(IntType n, std::set<IntType> &primeFactors) {
  DEBUG_FLAG(false);
  DEBUG("PrimeFactorize " << n);

  // primeFactors.clear();
  DEBUG("In PrimeFactorize n " << n);
  DEBUG("set size " << primeFactors.size());

  if (n == IntType(0) || n == IntType(1)) return;
  DEBUG("calling MillerRabinPrimalityTest(" << n << ")");
  if (MillerRabinPrimalityTest(n)) {
    DEBUG("Miller true");
    primeFactors.insert(n);
    return;
  }

  DEBUG("calling PrFact n " << n);
  IntType divisor(PollardRhoFactorization(n));

  DEBUG("calling PF " << divisor);
  PrimeFactorize(divisor, primeFactors);

  DEBUG("calling div " << divisor);
  // IntType reducedN = n.DividedBy(divisor);
  n /= divisor;

  DEBUG("calling PF reduced n " << n);
  PrimeFactorize(n, primeFactors);
}

// issue-881: make sure we don't overflow an IntType
template <typename IntType>
IntType FirstPrime(uint64_t nBits, uint64_t m) {
#if (NATIVEINT == 64 && !defined(HAVE_INT128)) || (NATIVEINT == 128)
  if ((typeid(IntType) == typeid(NativeInteger)) &&
      (nBits >= MAX_MODULUS_SIZE)) {
    PALISADE_THROW(math_error, "Requested modulus size " +
                                   std::to_string(nBits + 1) +
                                   " exceeds maximum allowed size " +
                                   std::to_string(MAX_MODULUS_SIZE));
  }
#endif
  try {
    DEBUG_FLAG(false);
    IntType r = IntType(2).ModExp(nBits, m);
    DEBUG("r " << r);
    IntType qNew = (IntType(1) << nBits);
    double d1 = qNew.ConvertToDouble();
    double d2 = std::pow(2, static_cast<double>(nBits));
    if ((qNew == IntType(0)) ||
        (std::llround(std::log2(d1)) < std::llround(std::log2(d2)))) {
      // too big for IntType
      PALISADE_THROW(math_error,
                     "FirstPrime parameters are too large for this integer "
                     "implementation");
    }
    IntType qNew2 = (r > IntType(0)) ? (qNew + (IntType(m) - r) + IntType(1)) :
                                       (qNew + IntType(1)); // if r == 0
    if (qNew2 < qNew) {
      PALISADE_THROW(
          math_error,
          "FirstPrime parameters overflow this integer implementation");
    }
    qNew = qNew2;

    while (!MillerRabinPrimalityTest(qNew)) {
      qNew2 = qNew + IntType(m);
      if (qNew2 < qNew) {
        PALISADE_THROW(math_error, "FirstPrime overflow growing candidate");
      }
      qNew = qNew2;
    }

    return qNew;
  } catch (...) {
    PALISADE_THROW(math_error, "FirstPrime math exception");
  }
  //  try {
  //    DEBUG_FLAG(false);
  //    IntType r = IntType(2).ModExp(nBits, m);
  //    DEBUG("r "<<r);
  //    IntType qNew = (IntType(1) << nBits);
  //    double d1 = qNew.ConvertToDouble();
  //    double d2 = std::pow(2, double(nBits));
  //    if ((qNew == IntType(0)) || (std::llround(std::log2(d1)) <
  // std::llround(std::log2(d2)))) {
  //      // too big for IntType
  //      PALISADE_THROW(math_error, "FirstPrime parameters are
  // too large for this integer implementation");
  //    }
  //    IntType mint = IntType(m);
  //    qNew -= mint + r - IntType(1);
  //
  //    while (!MillerRabinPrimalityTest(qNew)) {
  //      qNew -= mint;
  //    }
  //    return qNew;
  //  } catch (...) {
  //    PALISADE_THROW(math_error, "FirstPrime math exception");
  //  }
}

template <typename IntType>
IntType NextPrime(const IntType &q, uint64_t m) {
  IntType M(m);
  IntType qOld(q), qNew(q);

  do {
    qNew += M;
    if (qNew < qOld)
      PALISADE_THROW(math_error, "NextPrime overflow growing candidate");
  } while (!MillerRabinPrimalityTest(qNew));

  return qNew;

  return qNew;
}

template <typename IntType>
IntType PreviousPrime(const IntType &q, uint64_t m) {
  IntType M(m);
  IntType qNew = q - M;

  while (!MillerRabinPrimalityTest(qNew)) {
    qNew -= M;
  }

  return qNew;
}

template <typename IntType>
IntType NextPowerOfTwo(const IntType &n) {
  usint result = ceil(log2(n));
  return result;
}

/*Naive Loop to find coprimes to n*/
template <typename IntType>
std::vector<IntType> GetTotientList(const IntType &n) {
  std::vector<IntType> result;
  IntType one(1);
  for (IntType i = IntType(1); i < n; i = i + IntType(1)) {
    if (GreatestCommonDivisor(i, n) == one) result.push_back(i);
  }

  return result;
}

/* Calculate the remainder from polynomial division */
template <typename IntVector>
IntVector PolyMod(const IntVector &dividend, const IntVector &divisor,
                  const typename IntVector::Integer &modulus) {
  usint divisorLength = divisor.GetLength();
  usint dividendLength = dividend.GetLength();

  IntVector result(divisorLength - 1, modulus);
  usint runs = dividendLength - divisorLength + 1;  // no. of iterations

  // Precompute the Barrett mu parameter
  auto mu = modulus.ComputeMu();

  auto mat = [](const typename IntVector::Integer &x,
                const typename IntVector::Integer &y,
                const typename IntVector::Integer &z,
                const typename IntVector::Integer &mod,
                const typename IntVector::Integer &muBarrett) {
    typename IntVector::Integer result(z.ModSub(x * y, mod, muBarrett));
    return result;
  };

  IntVector runningDividend(dividend);

  usint divisorPtr;
  for (usint i = 0; i < runs; i++) {
    typename IntVector::Integer divConst(runningDividend.at(
        dividendLength - 1));  // get the highest degree coeff
    divisorPtr = divisorLength - 1;
    for (usint j = 0; j < dividendLength - i - 1; j++) {
      if (divisorPtr > j) {
        runningDividend.at(dividendLength - 1 - j) =
            mat(divisor.at(divisorPtr - 1 - j), divConst,
                runningDividend.at(dividendLength - 2 - j), modulus, mu);
      } else {
        runningDividend.at(dividendLength - 1 - j) =
            runningDividend.at(dividendLength - 2 - j);
      }
    }
  }

  for (usint i = 0, j = runs; i < divisorLength - 1; i++, j++) {
    result.at(i) = runningDividend.at(j);
  }

  return result;
}

template <typename IntVector>
IntVector PolynomialMultiplication(const IntVector &a, const IntVector &b) {
  usint degreeA = a.GetLength() - 1;
  usint degreeB = b.GetLength() - 1;

  usint degreeResultant = degreeA + degreeB;

  const auto &modulus = a.GetModulus();

  IntVector result(degreeResultant + 1, modulus);

  for (usint i = 0; i < a.GetLength(); i++) {
    for (usint j = 0; j < b.GetLength(); j++) {
      const auto &valResult = result.at(i + j);
      const auto &valMult = a.at(i) * b.at(j);
      result.at(i + j) = (valMult + valResult).Mod(modulus);
    }
  }

  return result;
}

template <typename IntVector>
IntVector GetCyclotomicPolynomial(usint m,
                                  const typename IntVector::Integer &modulus) {
  auto intCP = GetCyclotomicPolynomialRecursive(m);
  IntVector result(intCP.size(), modulus);
  for (usint i = 0; i < intCP.size(); i++) {
    auto val = intCP.at(i);
    if (intCP.at(i) > -1) {
      result.at(i) = typename IntVector::Integer(val);
    } else {
      val *= -1;
      result.at(i) = modulus - typename IntVector::Integer(val);
    }
  }

  return result;
}

template <typename IntVector>
typename IntVector::Integer SyntheticRemainder(
    const IntVector &dividend, const typename IntVector::Integer &a,
    const typename IntVector::Integer &modulus) {
  auto val = dividend.at(dividend.GetLength() - 1);

  // Precompute the Barrett mu parameter
  auto mu = modulus.ComputeMu();

  for (int i = dividend.GetLength() - 2; i > -1; i--) {
    val = dividend.at(i) + a * val;
    val = val.Mod(modulus, mu);
  }

  return val;
}

template <typename IntVector>
IntVector SyntheticPolyRemainder(const IntVector &dividend,
                                 const IntVector &aList,
                                 const typename IntVector::Integer &modulus) {
  IntVector result(aList.GetLength(), modulus);
  for (usint i = 0; i < aList.GetLength(); i++) {
    result.at(i) = SyntheticRemainder(dividend, aList.at(i), modulus);
  }

  return result;
}

template <typename IntVector>
IntVector PolynomialPower(const IntVector &input, usint power) {
  usint finalDegree = (input.GetLength() - 1) * power;
  IntVector finalPoly(finalDegree + 1, input.GetModulus());
  finalPoly.at(0) = input.at(0);
  for (usint i = 1; i < input.GetLength(); i++) {
    finalPoly.at(i * power) = input.at(i);
  }
  return finalPoly;
}

template <typename IntVector>
IntVector SyntheticPolynomialDivision(
    const IntVector &dividend, const typename IntVector::Integer &a,
    const typename IntVector::Integer &modulus) {
  usint n = dividend.GetLength() - 1;
  IntVector result(n, modulus);

  auto mu = modulus.ComputeMu();

  result.at(n - 1) = dividend.at(n);
  auto val(dividend.at(n));
  for (int i = n - 1; i > 0; i--) {
    val = val * a + dividend.at(i);
    val = val.Mod(modulus, mu);
    result.at(i - 1) = val;
  }

  return result;
}

}  // namespace lbcrypto
